library(haven)
library(ggplot2)
library(tidyr)
library(dplyr)
library(stringr)

setwd("C:/Users/jonves/Dropbox/Papers/SSP-Instability-Opinion/replication_data/")

growth <- function(x){x/lag(x)-1}

pwt <- read_dta("data/pwt90.dta")
pwt2gwno <- read.csv("data/pwt2gwno.csv", stringsAsFactors = FALSE, strip.white = TRUE)

pwt <- pwt[which(!is.na(pwt$rgdpna)),]

# Exclude ARE, QAT, KWT, BRN, SAU, BHR because they are extreme outliers, especially before, during and after the oil crisis in the 70s.
oil <- c("ARE", "QAT", "KWT", "BRN", "SAU", "BHR")


df <- group_by(pwt, countrycode) %>%
  filter(min(year)<=1970) %>%
  filter(!countrycode %in% oil) %>%
  filter(year >=1970) %>%
  mutate(grwt = growth(rgdpna/pop)) %>%
  summarise(mean_grwt = mean(grwt, na.rm=TRUE),
            gdpcap1970 = rgdpna[which(year==1970)]/pop[which(year==1970)],
            gdpcap2014 = rgdpna[which(year==2014)]/pop[which(year==2014)])



df$y <- df$mean_grwt
df$x <- df$gdpcap1970

p1 <- ggplot(df, aes(x=x, y=y)) + geom_point() +
  geom_smooth(data=df, aes(x=x, y=y), inherit.aes=FALSE, method = "lm", se=FALSE, color="red", formula = y ~ x) +
  ylab("Average per capita growth between 1970 and 2014") + xlab("GDP/capita in 1970") +
  scale_x_log10() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_bw() +
  theme(strip.text.x = element_text(size=12),
        axis.title = element_text(size=12),
        axis.text = element_text(size=12),
        legend.position = "bottom",
        legend.text = element_text(size=12),
        legend.title = element_text(size=12)) +
  guides(color=guide_legend(nrow=2,byrow=TRUE))


df$y <- log(df$gdpcap2014/df$gdpcap1970) # This produces the same plot (plots average growth rates instead of average growth).
ggplot(df, aes(x=x, y=y)) + geom_point() +
  geom_smooth(data=df, aes(x=x, y=y), inherit.aes=FALSE, method = "lm", se=FALSE, color="red", formula = y ~ x) +
  ylab("Growth rate between 1970 and 2014") + xlab("GDP/capita in 1970") +
  scale_x_log10() +
  theme_bw() +
  theme(strip.text.x = element_text(size=12),
        axis.title = element_text(size=12),
        axis.text = element_text(size=12),
        legend.position = "bottom",
        legend.text = element_text(size=12),
        legend.title = element_text(size=12)) +
  guides(color=guide_legend(nrow=2,byrow=TRUE))



ggsave("results/figure1a.tiff", plot=p1, width=6, height=6, dpi=300, compression="zip")